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Abstract 

A domain decomposition method is developed 
for solving the aerodynamic/ aeroacoustic problem 
of an airfoil in a vortical gust. The computa- 
tional domain is divided into inner and outer regions 
wherein the governing equations are cast in different 
forms suitable for accurate computations in each re- 
gion. Boundary conditions which ensure continuity 
of pressure and velocity are imposed along the inter- 
face separating the two regions. A numerical study is 
presented for reduced frequencies ranging from 0.1 
to 3.0. It is seen that the domain decomposition 
approach is far superior to the conventional single 
domain approach in providing robust and grid inde- 
pendent solutions. 

I. Introduction 

Many flow fields that occur in aerospace ap- 
plications involve upstream flow disturbances which 
propagate downstream, interact with structural 
components, and radiate sound. Typical examples 
include turbulent flow past a wing, unsteady flow 
around a propeller blade, and unsteady flow through 
a row of stator blades. 

The governing equations for such flows are the 
unsteady Navier-Stokes equations. However, viscous 
effects are often confined to small regions of the flow, 
and the unsteady Euler equations can be solved in- 
stead. If one further assumes that the converted 
disturbances are not too large, and that the flow 
moves at high speed, then one can invoke the “rapid 
distortion” 1,2 approximation and solve the linearized 
unsteady Euler equations. In this case, one obtains 
the zeroth-order steady mean flow first, and then ob- 
tains the unsteady flow as a first-order perturbation. 
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When the mean flow is both inviscid and irro- 
tational, the steady velocity can be expressed as the 
gradient of a potential, and the problem for deter- 
mining the unsteady flow can be reduced to solving a 
single convective wave equation with a dipole source 
term. This was first shown by Goldstein 3 , who de- 
composed the unsteady velocity into the sum of a 
potential component V</, and a vortical component 
u^\ so that u{x,t) = V(j) + u < ' 1 \ u ^ is a known 
function of the upstream velocity disturbances and 
the mean flow Lagrangian coordinates. (f> satisfies 


A / 1 A/ 1 
Dt c 0 2 Dt 


— V • (poV0) 
Po 


-V-( W a (/) ), 

Po 


where Co is the mean flow speed of sound, po is the 
mean flow density, and jx the convective deriva- 
tive based on the mean flow velocity. The unsteady 
pressure is given by p = —po~px- Goldstein’s for- 
mulation thus reduces the linearized Euler equations 
to a single convective wave equation. 

For most aerodynamic flows, there will be a 
frontal stagnation point or line where the mean ve- 
locity vanishes and the Lagrangian coordinates be- 
come singular. In this case, u D') also becomes sin- 
gular and remains so along the body surface. This 
makes it difficult to use Goldstein’s formulation di- 
rectly for numerical calculations, since the potential 
velocity component V<^ must cancel the singular be- 
havior of wA to satisfy the impermeablity condition. 

Atassi and Grzedzinski 4 proposed a modified 
splitting of the unsteady velocity which removes the 
singular and indeterminate character of the vorti- 
cal velocity on the body surface. Here the un- 
steady velocity is decomposed into the sum of an un- 
known potential component V<^, and a known vorti- 
cal component u^ R \ where has zero normal and 
streamwise velocity components on the body surface. 
The potential 4> satisfies Goldstein’s convective wave 
equation with a modified source term, 


A , 1 Dq<p 


— V • OoV(/) 
Po 


-V-(p 0 u (R) ). 

Po 
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The surface boundary condition for <f> is no longer 
singular and is in fact just V(j> ■ n = 0. 

In a series of papers 5-10 , Scott and Atassi pre- 
sented the first numerical implementation of Atassi 
and Grzedzinski’s linearized formulation for the case 
of unsteady vortical flow past a single airfoil. The 
main objective in this effort was to develop an aero- 
dynamic solver, the GUST3D code, which could 
accurately predict the airfoil unsteady response to 
three-dimensional, periodic vortical gusts. 

Another series of papers 11-15 focused on the 
far-held aeroacoustic response. It was noted in these 
papers that, although the GUST3D code produced 
an accurate near-held solution, there was a signifi- 
cant loss of accuracy in the far held. As a result, 
these papers concentrated on the development of 
Kirchoff methods to extend the GUST3D mid-held 
solution to the far held. Comparison of Kirchoff re- 
sults with analytical solutions showed this to be a 
promising approach. 

Building on this work, Scott 16 presented a se- 
ries of benchmark solutions for CAA code vali- 
dation at the recent Third Computational Aeroa- 
coustics Workshop on Benchmark Problems. Com- 
parison with results from nonlinear time-marching 
codes 16-21 showed good agreement on the airfoil sur- 
face and good agreement in the far held for low re- 
duced frequencies. However, there were some dis- 
crepancies in the high reduced frequency compar- 
isons. 

An extensive evaluation of the GUST3D code 
indicated that, as the reduced frequency increases, 
the source term calculated by the code grows rapidly 
in the far held and oscillates, making it difficult to 
obtain an accurate solution. 

The purpose of this paper is to present a new 
domain decomposition approach which largely cor- 
rects this problem. The basic idea is to divide the 
how held into inner and outer regions, and to use the 
Atassi-Grzedzinski formulation in the inner region 
and the Goldstein formulation in the outer region. 
It will be seen that this approach leads to substantial 
improvements in accuracy, both on the airfoil and in 
the far held. 


(Aq, /c2, £3) is the wave number vector, and a and k 
satisfy a ■ k — 0 to ensure that the continuity equa- 
tion is satishecl. 

Let the how held be represented by 


U(x, t) = Up(x) + u(x, t) 

(2.2) 

p{x,t) = p 0 (x) + p'(x,t) 

(2.3) 

p{x,t) = po(x) + p'{x,t) 

(2.4) 

s(x, t) = So + s'(x, t ) 

(2.5) 


where the entropy So is constant, and u, p\ p', and 
s' are the unsteady perturbation velocity, pressure, 
density and entropy, respectively. Quantities with 
“0” subscripts are the steady mean how quantities 
which are assumed to be known. 

Substituting (2.2) - (2.5) into the nonlinear Eu- 
ler equations and neglecting products of small quan- 
tities, one obtains the linearized continuity, momen- 
tum, and entropy conservation equations 

+ p'V • U 0 + V • (p 0 u) = 0 (2.6) 


Po( 


D 0 u 

Dt 


+ u ■ V/T 0 ) + p'U Q • W 0 = -Vp' 


Dps ' 
Dt 


= 0, 


(2.7) 

( 2 . 8 ) 


where ^ f + Up • V. 

Let the how held be divided into inner and outer 
regions, as shown in Figure 1. In the outer region, 
let the unsteady velocity be decomposed according 
to Goldstein’s velocity splitting, 


u(x,t) = V(j) 0 + u( G \ (2.9) 


where we define = u and where the “o” 

subscript denotes the outer region. Equations (2.6) 
- (2.8) are then reduced to 


C<j>o = -V-(p 0 w (G) ) (2.10) 

Pp 


II. Mathematical Formulation 


Governing Equations 

Consider an airfoil with chord length c in a how 
with uniform upstream velocity /Too in the x\ direc- 
tion. Let the fluid be an ideal gas which is inviscid 
and non-heat-conducting. Far upstream, let 


W = ae ik < s ~ iU °°V 


denote a small amplitude gust, where i is a unit 
vector in the x\ direction. Here a = (01,02,03), 
where the amplitude |a| satisfies |a| << /Too, k = 


r _ A), 1 D 0 1 -• 

£ — TTA - 2777) “ — V - poV . 
Dt co Dt po 

The unsteady pressure is given by 


P = (2-12) 

In the inner region, let the velocity be decom- 
posed according to the Atassi-Grzedzinski velocity 
splitting, 


T(x, t) = + ' 
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where the subscript denotes the inner region. 
Equations (2.6) - (2.8) then reduce to 

£</>/ = — V- (p 0 n (ii) ). (2.14) 

Po 

For flows with no upstream entropy disturbances, 
u and u( G ' ) are related by 3 4 

u^ R) = n (G) + V4>, (2.15) 

where <f> is a function which is constructed to cancel 
the singularity in u^ G ' ) on the airfoil surface, so that 
(j), has a regular boundary condition. It should be 
noted that (j) has no pressure associated with it, so 
that the unsteady pressure in the inner region is just 


Boundary Conditions 

At the airfoil surface, the normal velocity com- 
ponent must vanish. Since u^ R> has zero normal and 
streamwise velocity components along the airfoil sur- 
face, the airfoil boundary condition is just 

V(t>j-n = 0. (2.17) 

Across the wake, in both the inner and outer 
regions, the pressure and normal velocity must be 
continuous. For the pressure, this leads to 


where the spatial derivatives in (2.21) are under- 
stood to be one-sided. 

For the velocity, one uses (2.9) and (2.13) to- 
gether with (2.15) to obtain 

V0o = V0J + V0. (2.22) 

Taking the dot product of each term in (2.22) with 
the interface unit normal ft, one obtains the conti- 
nuity of normal velocity condition 

V<j) 0 -n = V</> 7 • n + V(f>-n. (2.23) 

Alternatively, one can integrate (2.22) to obtain the 
jump condition 

<t>o ~ <t> i = 4>- (2-24) 

The pressure and normal velocity conditions, 
equations (2.21) and (2.23), represent a consistent 
set of conditions for non-overlapping domain de- 
composition. Equation (2.23) is a Neumann con- 
dition, while equation (2.21) is a linear combination 
of Dirichlet and Neumann conditions. They always 
form a linearly independent and therefore robust set 
of interface conditions for domain decomposition. 

On the other hand, at an interface location 
where the mean velocity becomes tangent to the in- 
terface, equations (2.21) and (2.24) may not be lin- 
early independent. This could lead to difficulties for 
some interface configurations. 


P+ - P- = 0 ) = 0, (2.18) 

where “+” and subscripts denote quantities 
above and below the wake, respectively. (The inner 
and outer subscripts are omitted here for simplicity.) 
The inner and outer potentials must then satisfy 


A>(A0) 

Dt 


(2.19) 


where A (j) denotes the jump in <f> across the wake. 
For continuity of normal velocity, the potentials 
must satisfy 


V0 + • ft = V<^_ • ft. (2.20) 

On the outer grid boundary, a radiation bound- 
ary condition must be imposed. This will be dis- 
cussed in the next section. 

Finally, interface conditions must be specified at 
the boundary separating the inner and outer regions. 
Analogous to the wake, the pressure and velocity 
should be continuous. For the pressure, this requires 
that 

Doto = Do fa ( 2 . 21 ) 

Dt Dt K ’ 


III. Numerical Implementation 

For the most general case, in which the steady 
velocity Uq(x) varies spatially as a function of x, 
the right hand sides of (2.10) and (2.14) are ex- 
pressed as functions of the mean flow Lagrangian 
coordinates. One therefore introduces the variables 
(Xi,X 2 , X 3 ) = X, where 


A 2 


T 0 

Poo b'oo 


and 


A3 — X3, 


(3.1) 

(3.2) 


where To is the stream function of the mean flow and 
X 3 is the spatial coordinate in the spanwise direction. 
The component X\ is given by 


Ai = U 0 0 A, (3.3) 


where A is the Darwin-Lighthill “drift” 
function 22,23 , which can be expressed in terms of 4>o 
and Tq as 


A = 


To 

t/00 



Tj - 2 ) d ®0- 

U OO 


(3.4) 
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The integration is carried out on To = constant (i.e., 
a fixed streamline). 

For flows with no upstream entropy fluctua- 
tions, Goldstein’s vortical velocity may then be ex- 
pressed as 

u (G) = [V(o • (3.5) 

and the Atassi-Grzedzinski vortical velocity is 

u {R) = [V(a • x)]e* (X '" rc/ “ t) + V0. (3.6) 


Equations (2.10) and (3.11) are most conve- 
niently solved in the frequency domain using the 
(T 0 ,To) orthogonal curvilinear coordinate system, 
where To and To are the mean flow potential and 
stream functions. 

We now assume that all variables have been 
nondimensionalized as in [7,10]. The normalized 
wave number k± denotes the reduced frequency, and 
the free stream Mach number is denoted by 
We transform into the frequency domain in the in- 
ner and outer regions using 


The function 0 is constructed to cancel the singular- 
ity in tS G ' > on the airfoil surface, and is given by 4 


0 = 



where 


a 2 k\ — a\k 2 1 — e lfe2A ’ 2 
1 + iaoUooki k 2 


faik^X-iU^t) 

(3.7) 


( dUo. 

ao = 


-1 


(3.8) 


Here n denotes the direction of the outward unit 
normal, and S denotes the stagnation point near the 
airfoil leading edge. 

Now from (3.5) and (2.1), one can show that as 
X\ — > —oo, — > Uqo. Since we must also have 

u(x, t) — > Uoo at upstream infinity, it follows that 
in the outer region 0 O must satisfy V0 O — > 0 as 
X\ — > —oo. This ensures that <j> 0 has outgoing wave 
behavior at infinity. 

On the other hand, in the inner region, it fol- 
lows from (2.1), (2.13) and (3.6) that (f), must satisfy 
V0j — > — V0 as one moves toward upstream infin- 
ity. It is necessary, therefore, to replace 0* with a 
function whose gradient vanishes as x\ — > — oo. This 
will ensure that the new potential has outgoing wave 
behavior, and reduce any incompatibility across the 
interface separating the inner and outer regions. 

Following the formulation presented in [7,10], 
we introduce the potential functions 0 i and 02 , 
where 

<t>i = <t > l - 02, (3.9) 


and 0 2 is a known function which is constructed such 
that 


[02 — 0 | — > 9 as r — > oo, (3.10) 

where r denotes polar distance. This ensures that 
V0i — > 0 in the far field. 

Upon substituting (3.9) into (2.14), the inner 
governing equation becomes 

£0! = -v-(p 0 u (fl) ) + £ 02 , (3.11) 

Po 


01 = ( p I e -ikit+ik 3 x 3 (3.12) 

and 

0o = v 0 e- iklt+ik3X3 . (3.13) 

Transformation into computational coordinates 
is then accomplished as follows. First, introduce 
Prandtl-Glauert coordinates (T,T) by 

T = T 0 (3.14a) 

T = /3ooT 0 , (3.14b) 


where (3 Then introduce new depen- 
dent variables 0j and ip 0 , where 


Pi = ipie 

(3.15) 

Po = 0o e ~ lK ° $ 

(3.16) 

h 

Ao= 0 2 • 

(3.17) 


Finally, transform T and T into computational co- 
ordinates using 

T = a* cos(7r?7) cosh(7r£) 

(3.18a) 

T = a* sin(7rr7) sinh(7r£), 

(3.18b) 


where a* is a known constant. The inner governing 
equation then becomes 



<9 2, 0j 

~w 




+ A t J M ^ +Tl ^ + T 2 ^ + 


n 9 2 0 1 


(9 2 0, d 2 0 f 

4 c>77 2 5 dr]d^ 


= s„ 


(3.19) 


where the right hand side consists of known func- 
tions. The governing equation in the outer region 
remains that given in (2.10) 


where J, Ai, and T) ... T 5 are known functions, and 
S, is the source term. Similarly, the outer governing 
equation is 
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IV. Numerical Results 


-Pi 


r a^o 

. dr/ 2 


d 2 j> 0 
d £ 2 




+ Ai j ( 77 , £) %i) 0 + Ti + 

aVo . r 

4 dr] 2 5 dqdi 


■ r dip 0 

T2 ^f 

= S 0 . 


+ T 3 


d 2 ^o 

d? 


(3.20) 


Equations (3.19) and (3.20) are implemented 
using nine-point central differencing which is second- 
order accurate. Each equation is imposed at interior 
grid points within its respective region. 

At the interface separating the two regions, one 
must ensure that the pressure and velocity are con- 
tinuous, as discussed in the previous section. This 
is accomplished by using a row of coincident grid 
points along the interface, with one set of points 
belonging to the inner region and the other set be- 
longing to the outer region. (See Figures 2a and b.) 
Pressure continuity is satisfied by imposing condi- 
tion (2.21). Velocity continuity is satisfied by impos- 
ing either (2.23) or (2.24). In calculations to date, 
conditions (2.23) and (2.24) have been found to give 
nearly identical results. 

For wake grid points, continuity of pressure and 
normal velocity are enforced by way of equations 
(2.19) and (2.20). Equation (2.19) is imposed in 
integral form for every wake point on the upper side. 
Equation (2.20) is imposed using three-point, one- 
sided differencing for every wake point on the lower 
side. 

On the outer grid boundary, we impose the 
Bayliss-Turkel 24 radiation boundary condition of or- 
der 1. This condition is applied to the unsteady 
pressure, and can be written 


(£ - b )(m- a ) = »• < 3 - 2I > 

where 

A = % (3-22) 

and 

B = ik - (3.23) 

2 r 

Here k is the Helmholtz constant which is defined by 


k 2 = 


/fclMoo\2 /l- 3 \ 2 

l /& ) V/3o J ' 


(3.24) 


Condition (3.21) has proven to be both accurate and 
computationally efficient 25 . 


In this section, we compare numerical results 
using the new domain decomposition approach ver- 
sus the original single domain approach. All cal- 
culations are for a 12% thick, symmetric Joukowski 
airfoil in a 2-D gust propagating at 45°, i.e., = k\. 

The airfoil has zero degrees angle of attack and no 
mean loading. The Mach number is 0.5. We consider 
reduced frequency values k\ = 0.1, 1.0, 2.0, and 3.0 
(with normalization based on the half chord). The 
gust amplitude is taken to be 2% of the free stream 
velocity. 

In the results that follow, we present RMS pres- 
sure on the airfoil surface and acoustic intensity in 
the far field. For each reduced frequency, we present 
results from a series of calculations on five different 
grids. Each grid differs only in the location of its 
outer grid boundary. The mesh spacing is the same 
for all five grids, with uniform 77 spacing and vari- 
able £ spacing. The £ spacing provides 24 points per 
gust wavelength. Our main objective is to assess the 
ability of each approach to give a consistent solution 
which does not depend on the outer grid boundary 
location. 

Figures 3 and 4 show the RMS pressure on the 
airfoil surface for the low frequency case k\ = 0.1. 
The legend at the top of the figure indicates the dis- 
tance (in gust wavelengths) to the outer grid bound- 
ary for each grid. This is specified in terms of the 
GUST3D parameter “nwaves”. The results in these 
figures show that the airfoil pressure is indeed grid 
independent for each approach. 

Figures 5 and 6 show the corresponding acous- 
tic intensity on a circle of radius two chord lengths, 
centered about the airfoil center. (This circle lies 
within the inner region for all results presented in 
this paper. In general, there is no relationship be- 
tween the location of the circle and the location of 
the interface.) Each figure shows some sensitivity 
of the far-field pressure to the location of the outer 
grid boundary, with the domain decomposition re- 
sults being slightly less sensitive. 

We should note that the values of nwaves that 
were used in our calculations were designed to op- 
timize accuracy for each reduced frequency and for 
each computational approach. All figures show re- 
sults for five consecutive values of nwaves, where 
nwaves was incremented by 0.5. Each figure shows 
the best set of five consecutive results available for 
that case. This is the reason for the different values 
of nwaves that are shown. 

In Figures 7-10, we present results for the mid- 
frequency case, k\ = 1.0. Here the single domain 
results begin to show significant sensitivity to the 
change in grid, even on the airfoil surface. On the 
other hand, the domain decomposition results are 
very nearly grid independent on the airfoil and ac- 
ceptably grid independent in the far field. 
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In Figures 11 - 18, we present results for the 
relatively high reduced frequencies of k\ = 2.0 and 
k\ = 3.0. At these higher frequencies, the single do- 
main results deteriorate markedly due to the rapidly 
growing and oscillating source term, as shown in Fig- 
ure 19. The domain decomposition results, however, 
give an acceptably grid independent solution both on 
the airfoil and in the far field. 

We should point out, however, that the domain 
decomposition results are somewhat sensitive to the 
location of the interface separating the inner and 
outer regions, especially at high frequencies. We 
should note in addition that grid independence by 
itself does not imply accuracy. For this reason, the 
k\ = 2.0 and k\ = 3.0 results should be considered 
preliminary at this time. 

Summary 

In this paper we have presented a new domain 
decomposition approach for the single airfoil gust re- 
sponse problem. We divide the flow field into inner 
and outer regions, and use the Atassi-Grzedzinski 
linearized Euler formulation in the inner region, and 
Goldstein’s linearized Euler formulation in the outer 
region. This approach uses each formulation where 
it is best suited. In the inner region, the Atassi- 
Grzedzinski formulation cancels the singularity in 
Goldstein’s vortical velocity, and provides a bound- 
ary value problem with regular boundary conditions. 
In the outer region, far away from the airfoil singu- 
larity, Goldstein’s formulation provides a boundary 
value problem which is better suited for wave prop- 
agation in an open domain. Numerical results show 
that the single domain approach is very sensitive to 
the location of the outer grid boundary, and is un- 
able to provide a consistent or grid independent so- 
lution at the higher reduced frequencies. On the 
other hand, the domain decomposition approach is 
largely insensitive to the location of the outer grid 
boundary, and provides an acceptably grid indepen- 
dent solution for reduced frequencies ranging from 

0.1 to 3.0. 
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